Scaling Properties of Flexible Membranes from Atomistic Simulations: 

Application to Graphene 
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Structure and thermodynamics of crystalline membranes are characterized by the long wavelength 
behavior of the normal-normal correlation function G(q). We calculate G(q) by Monte Carlo and 
Molecular Dynamics simulations for a quasi-harmonic model potential and for a realistic potential 
for graphene. To access the long wavelength limit for finite-size systems (up to 40000 atoms) we 
introduce a Monte Carlo sampling based on collective atomic moves (wave moves). We find a 
power-law behaviour G(q) oc q~ 2+v with the same exponent r\ « 0.85 for both potentials. This 
finding supports, from the microscopic side, the adequacy of the scaling theory of membranes in the 
continuum medium approach, even for an extremely rigid material like graphene. 

PACS numbers: 63.20.Ry, 68.60.Dv, 81.05.Uw, 05.10.Ln 



Collective phenomena involving infinitely many de- 
grees of freedom are often characterized by scaling laws 
with power-law behavior of correlation functions. In 
three dimensional systems, this behavior occurs only at 
critical points [l], 0, ||. I n two dimensions (2D) the 
situation is different, and a whole temperature interval 
with "almost broken symmetry" and power-law decay of 
correlation functions frequently appears, the Kosterlitz- 
Thouless (KT) transition in 2D superfluids and super- 
conductors [H being a prototype example. Existence of 
real long range order, where correlation functions remain 
non-zero in the limit of infinite distance, is forbidden in 
such cases by the Mermin- Wagner theorem [s| due to 
the divergence of the contribution of soft modes to rel- 
evant thermodynamic properties. The theory of flexible 
membranes embedded in higher dimensions is an im- 
portant part of the statistical mechanics of 2D systems. 
Here, we investigate the scaling behavior of crystalline 
flexible membranes by means of atomistic simulations, 
using graphene 0, BCa], the simplest known membrane, 
as an example. 

In the flat phase, the membrane in-plane and out-of- 
plane displacements are parametrized by a D-component 
'stretching' phonon field u a (x), a = 1...D, and by a 
d c = d — D component out-of-plane height fluctuation 
ft(x), where d is the space dimension and D is the mem- 
brane dimension. Softening of bending modes makes this 
situation very similar to the KT model. A minimal phc- 
nomenological model for membranes isiust the elasticity 
theory described by the Hamiltonian p, [l(| : 

H=\Jd D x (n(V 2 h) 2 + pulp + ^L) (1) 

where k, fi and A are bending rigidity, shear modulus and 
Lame coefficient and 
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is the strain tensor. In harmonic approximation, by ne- 
glecting the last, non-linear, term in Eq. (JSJ, the bending 
(h) and stretching (u) modes are decoupled. 

The Hamiltonian ([I]) is quadratic in the phonon de- 
grees of freedom u which can be eliminated by Gaus- 
sian integration [|| [To| . In this way, the Hamiltonian 
can be rewritten only in terms of the Fourier compo- 
nents of the height h as the sum of a harmonic bending 
energy, quadratic in h, and an anharmonic energy, quar- 
tic in h, that results from the coupling of bending and 
stretching modes • If one neglects the latter term, the 
membrane becomes crumpled at any finite temperature 
with, for D = 2, the mean square height fluctuations 
(h 2 ) ~ l? and normal- normal correlation functions that 
diverge logarithmically at large distances. Nelson and 
Peliti [10| suggested that the above anharmonic term 
stabilizes the flat phase at least at temperatures much 
smaller than k. This flat phase is described by an ef- 
fective bending rigidity K,(q) ~ q~ ri and effective elas- 
tic moduli with power-law dependencies on q that par- 
tially suppress long wavelength bending fluctuations. As 
a result, the normal-normal correlation function remains 
finite, although (h 2 ) still diverges as (h 2 ) ~ L 2C > with 
C = 1 — 77/2 Thus, the flat phase is not truly flat, 
but still exhibits rather strong corrugation. The model 
(fTj), which is called the model of phantom membranes, 
has a transition to a crumpled phase at temperature of 
the order of k. The long wavelength limit was solved 
within the Self Consistent Screening Approximation in 
Ref. [ll] yielding r\ = 0.821. The discretized version of 
this model was investigated by Bowick et al. by means 
of Monte Carlo simulations giving 77 w 0.72 [12]. The 
term 'phantom' means that the model does not include 
self- avoidance, the natural condition of true physical sys- 
tems. It is assumed that self-avoidance removes the phase 
transition to the high temperature crumpled phase while 
the scaling properties of the 'flat' phase remain the same 
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as in phantom membranes. However, any kind of accu- 
rate statement about the model can be justified only in 
the limit d c — > oo and, strictly speaking, nothing can be 
said rigorously for the real case of d = 3, D = 2 and 
d c = l. 

To characterize the long wavelength limit of the height 
fluctuations we compare the results of atomistic simu- 
lations to the predictions of this theory for the normal- 
normal correlation functions G(q) = { | n g | 2 ) . Starting 
from Eq. (TTJ) an expression for G(q) has been given from 
general scaling considerations [(I [ToL [l5| in the form of an 
effective Dyson equation 



G- 1 (q) = G Q 1 (q) + S (q) 



(3) 



where Go is the value derived in harmonic approximation 

TN 



G (q) 



nS q 2 



and the self energy is 



ASq 2 ( q \ v 



(5) 



with N the number of atoms, So = L x L y /N the area 
per atom, T the temperature in units of energy, qg = 
2tt\/bJk, B the two-dimensional bulk modulus [13| and 
A an unknown numerical factor. 

Until recently, this phenomenological continuum model 
was the only way to describe the statistical mechanics 
of membranes since all known real membranes were 
too complicated for atomistic models. The situation has 
been changed drastically by the discovery of graphene Q 
which is the first example of a truly two-dimensional sys- 
tem (just one atom thick) and, thus, a prototype crys- 
talline membrane @, Q . The experimental observation of 
ripples in freely hanged graphene I14| stimulated a large 
theoretical activity [ll, M, lH Hill El HH . In partic- 
ular, using the accurate bond order potential for carbon 



LCBOPII 22j, we were able to simulate structural and 



thermodynamical properties of graphene at finite tem- 
peratures [Til . l2lj ] by straightforward Monte Carlo (MC) 
simulations. The simulations confirmed the existence 
of thermally induced intrinsic ripples at finite tempera- 
tures resulting in strong anharmonic effects. However, 
we found that the normal-normal correlation function 
could not be described by Eq. ^ over the whole range 
of q [Hj]. In fact, G(q) followed the power law result- 
ing from the harmonic approximation (phonon picture, 
r\ = 0) at large enough q, but, after bending, at smaller 
q's we found a drop of the correlation functions not com- 
patible with a power law. Our conjecture at that time 
was that the extreme rigidity of graphene could be the 
reason why it could not be described by the phenomeno- 
logical theory of membranes in a continuum medium ap- 
proach jiq |. However, we felt that this point deserved 
further investigation. Here, we focus on the low-g re- 
gion in order to establish firmly whether a scaling law 
exists and, if so, to determine the scaling exponent. To 



this purpose, we simulate large systems, introduce new 
MC moves for phase space sampling and examine more 
than one model of the interatomic forces, including a 
simple quasi-harmonic (QH) model that yields a not too 
rigid membrane and the extremely rigid case of graphene 
which is well described by LCBOPII. In addition, for the 
QH model we verified ergodicity of our MC simulations 
by comparing with Molecular Dynamics (MD) results. 

We begin by considering a relatively simple QH model 
with energy given by: 



U 



\ ( K r{r%j - r eq ) 2 + Kg Y (Vtjk - Veqf J 



where the summations over j and k are over the nearest 
neighbors of atom i, — cos 9ijk and y eq = cos" 



eg j 

(4) with r eq = 1.42 A and 9 eq — 2tt/3 the ground state 
equilibrium nearest neighbor distance and bond angle 
in graphene. The stretching and the bending force con- 
stants, K r = 22 eV A~ 2 and Kg = 4 eV, respectively, 
were chosen to yield elastic moduli for isotropic and uni- 
axial compressions equal to those for the LCBOPII [2l| . 

In Fig. Q] we show the function G(q)/N (dotted line) 
calculated by extensive standard Monte Carlo simula- 
tions in the canonical ensemble at 300 K for a system 
with N = 37888, L x = 314.82 A and L y = 315.24 A and 
periodic boundary conditions in the a;y-plane. Starting 
from the Bragg peak at q = 47r/(3r eg ) = 2.94 A" 1 and 
going towards lower q we find, first, the power law rj = 
due to the harmonic contribution, then, a smaller slope 
followed by a drop at the smallest q < 0.08 A -1 which 
corresponds to a wavelength of about 75 A. This drop 
is similar to the one mentioned above and found previ- 
ously in Ref. [TH with the LCBOPII for graphene. These 
results are obtained by averaging over many configura- 
tions in the canonical ensemble obtained by the ordinary 
MC procedure which is based on random displacements 
of randomly chosen individual atoms and volume (area) 
fluctuations with a Metropolis acceptance rule. By us- 
ing Eq. (|U) we find that the bending constant for the 
QH potential is k — 0.4 eV, much softer than the 1.1 eV 
appropriate for graphene [15|, due to neglected interac- 
tions beyond first neighbors. The observation that also 
the simple QH model shows a suppression of long wave- 
length excitations made us think of the possibility that 
standard MC is not an efficient sampling technique in 
this case. To resolve this issue we (i) extended our MC 
phase space sampling with a new type of collective trial 
events that we call 'wave moves' described below, and (ii) 
performed MD simulations for the QH model [23| , allow- 
ing a direct comparison with the MC results, with and 
without wave moves. The equivalence of time averages in 
MD simulations with ensemble averages in MC simula- 
tions guarantees that the system is in thermodynamical 
equilibrium (ergodic) . 

In Fig. [T] we compare the results of standard MC with 
the results obtained by MD and by MC with the addition 
of wave moves. The MD results coincide with the stan- 
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FIG. 1: (color online) Normal-normal correlation functions 
G(q)/N calculated for a graphene system with N = 37888 by 
ordinary MC simulations (red-dashed line), MD simulations 
and MC simulations with wave moves with the QH poten- 
tial. The dashed lines show the asymptotic harmonic behav- 
ior with power laws q~ 2 for large q and the long-wavelength 
limit q-' 2 -"' with 77 = 0.85. 
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FIG. 2: (color online) Normal-normal correlation functions 
G(q)/N calculated for three systems with N = 12096 (L x — 
177.08 A, L y = 178.92 A), N = 19504 (L x = 226.27 A, L y = 
225.78 A), N = 37888 (L x = 314.82 A, L y = 315.24 A) by 
MC simulations with wave moves with LCBOPII. The dashed 
lines show the asymptotic harmonic behavior with power laws 
q~ 2 for large q and the long- wavelength limit g - ' 2- ^ with 
77 = 0.85. The dashed-dotted line is G a of Eq. © with the 
coefficients fixed by the asymptotic behavior. One can see 
that the crossover is much sharper in the simulations. 



dard MC in the range where the latter is described by a 
power law, but does not show the drop at small q and 
keeps the same slope till the smallest possible q allowed 
by our finite size system. The results of MC simulations 
with wave moves coincide for all q's with those obtained 
by MD, implying that the system is in thermodynamic 
equilibrium. Both curves display a power-law behavior 
for the whole range of q in the long wavelength limit. A 
best fit of the data yields an exponent 77 = 0.85. 

A wave move consists of a transversal, wavelike dis- 
placement of all atoms in the z-direction, perpendicular 
to the graphene plane. For a given wavevector q there 
are two possible, linearly independent wave excitations, 



yielding z-coordinate displacements for all atoms i 

Azj = (0.5 — R)As : q cos(qrj) and 
Azi = (0.5 - R)A s , q sin(cp-i) 

where r; is the 3D position of atom i and R is a random 
number between and 1. The amplitude As, q is chosen 
such that the acceptance rate for such a wave move is 
between 0.4 and 0.5. The appropriate value of As lQ de- 
pends on the size of the 2D box S = L x L y and on the 
wavevector q (see below). 

Due to the periodic boundary conditions in the x- and 
y-directions the candidate wavevectors for wave moves 
can be restricted to a set on a 2D grid: 




with integer m x and m y . This set was further bounded 
by applying only wave moves of long wavelengths since 
short wavelengths are already efficiently sampled by the 
individual atom displacement trials. Hence, we consider 
a finite set of {m x , 77ij,)-pairs corresponding to q- vectors 
within a circular region with radius q max around q = 0. 
This set was kept constant during the entire simulation. 
We choose q ma x equal to the g-value below which G(q) 
starts to bend down in standard MC simulation. More 
precisely, we took q m ax — 0.16, corresponding to a mini- 
mal wavelength of 40 A. Since transversal phonon modes 
have quadratic dispersion uj{q) ~ q 2 , the energy change 
associated with a wave move behaves as AE wm ~ A 2 S q q 2 . 
Therefore, we took A$ } q = A$/q to obtain similar accep- 
tance rates for each of the allowed g-vectors, as was in- 
deed confirmed by our simulations. This choice leaves one 
adjustable parameter, ^5. For different system sizes, the 
appropriate As roughly scales as A/S, but a correction 
is required to fine-tune the acceptance rate. On average, 
a wave move was attempted every MC step by choosing 
randomly one of the 2N q possible waves. Here, N q is 
the number of allowed q- vectors (or (m x , m y )-pairs) and 
the factor 2 comes from the fact that each wavevector 
yields two possible waves: a sine and a cosine wave. An- 
other random number R € (0,1) was then pulled to fix 
the amplitude (0.5 — R)As/q. Following the Metropolis 
procedure, a wave move is always accepted if the energy 
change AE wm is negative, whereas for AE wm > it is 
accepted with probability P = exp(— (3AE wm ) 1 requiring 
another random number R' G (0, 1) to decide for accep- 
tance when R' < P or rejection when R' > P. 

MD simulations are much more demanding than MC 
simulations and are not within reach for the rather com- 
plex LCBOPII potential for the present system size. The 
previous results with the QH harmonic potential, how- 
ever, show that equilibrium can be reached using MC 
with wave moves. The correlation function G(q) calcu- 
lated by MC with wave moves for LCBOPII are shown in 
Fig. [U for three system sizes. Again, we see the crossover 
from the harmonic behavior to a power law with r; = 0.85 
up to the smallest wavevectors. The main difference with 
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FIG. 3: (a) Quadratic out-of-plane displacement (h 2 ), aver- 
aged over all particles, as a function of the MC step for the 
same systems as in Fig. [2] The dashed horizontal lines de- 
note ({h 2 )}, the average (h 2 ) over all MC steps except the 
first 3 x 10 s steps of equilibration, (b) {(h 2 )) as a function of 
the average linear system size L = ^JL X L V compared to the 
scaling law {(h 2 )) = CL 2 ~" with C=0. 00232 and 77 = 0.85 
(dashed line). Both axis are in logarithmic scale. 



the results obtained with the QH potential is that, due 
to a higher bending rigidity, the crossover between the 
two power laws is shifted to lower q values. Moreover, 
we note that for q > 1 A -1 there is a deviation from a 
power law behavior just before the Bragg peak. 

Finally, in Fig. [3Ja) we show the average out-of-plane 
displacement (h 2 ) corresponding to the simulations for 



LCBOPII of Fig. [2] which shows large fluctuations. In 
Fig. [3Kb) we plot the values of (h 2 ) averaged over all 
MC steps as a function of the system size in comparison 
with the expected scaling law exp(2 — 77). Although it 
would have been impossible to deduce the scaling expo- 
nent from the three points in Fig. E^b) due to the large 
error originating from the large fluctuations, these re- 
sults are certainly compatible with the scaling exponent 
r} found by a fit of G(q). With (h 2 ) = 1.65 A 2 for L = 315 
A and 77 = 0.85 we estimate \J (h 2 ) « 9 A for L=ljum, 
well in the range of measured values [l4[ • 

In summary, we have shown by atomistic simula- 
tions that, in thermodynamic equilibrium, crystalline 
membranes display a power-law scaling behavior of the 
normal-normal correlation function, in qualitative agree- 
ment with continuum medium theory. For different mod- 
els of interactions with different rigidities, we found the 
same exponent of anomalous bending rigidity r\ w 0.85. 
We have demonstrated that the efficiency of MC simu- 
lations for this type of systems can be greatly improved 
by introducing collective wave moves. On the basis of 
our results, we conclude that despite its extreme rigidity, 
graphene behaves as a prototype membrane opening new 
ways to study the intriguing physics of membranes on a 
system with well known interatomic interactions. 

This work is part of the research program of the 'Sticht- 
ing voor Fundamenteel Onderzoek der Materie (FOM)', 
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